# Replication Code for PGMs Paper


#---Packages---#
library(dplyr)
library(tidyverse)
library(ggplot2)
library(miceadds)
library(stargazer)
library(sandwich)
library(lmtest)
library(jtools)
library(sjPlot)
library(ggeffects)
library(fixest)
library(pglm)
library(plm)
library(margins)
library(broom)
library(estimatr)
library(huxtable)
library(ggpubr)
library(grid)


#---Data---#
data <- read.csv('full data all variables.csv') %>%
  select(-X) 

pgmcount <- read.csv('pgmdv.csv')

data <- data %>%
  left_join(pgmcount)


#--- Creating Table 1: Summary Statistics ---#

#-- read in data
table1 <- data %>%
  select(pgm2,log_pgm_count_stata,lag.roll_avg_hostlev,lag.roll_internal_conflict_intensity,log.lag.gdppc_i,lag.polity2,lagdemoc,lag.treaties_signed)

#-- create table
stargazer(table1, out="papersummarystats.tex",
          covariate.labels=c("Smart Bomb Acquisition (Binary)","Logged Count of Smart Bomb Acquisition","Lagged Average MID Hostility Level","Lagged UCDP Internal Conflict Intensity","Lagged Log GDP Per Capita","Lagged Polity2 Score","Lagged Democracy (Binary)", "Lagged ICRC Treaties Signed"),
          style = "ajps", summary.stat = c("n", "mean", "sd","min","max")
)


#--- Creating Figure 2 ---#


#----- Internal Threats Subfigure -----#

m1.2 <- glm(formula=pgm2~lag.roll_internal_conflict_intensity+lag.roll_avg_hostlev+lag.treaties_signed+lag.polity2+log.lag.gdppc_i #+t+t2+t3
            , data=data, 
            family=binomial)

summary(m1.2)


#-- getting the predicted values
internal <- ggpredict(m1.2, terms = "lag.roll_internal_conflict_intensity")

#-- plotting

tiff("Fig2Internal.tiff", units="in", width=5, height=5, res=300)

ggplot(internal, aes(x, predicted)) + 
  geom_line(color="#259797") +
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high), fill="#259797", alpha=0.15) +
  theme_bw() +
  ggtitle("Internal Threats") +
  theme(plot.title = element_text(hjust = 0.5, family="Times"), axis.title = element_text(family="Times"), axis.text = element_text(family="Times")) +
  labs(x= "UCDP Internal Conflict Intensity", y = "Predicted Probability of Smart Bomb Acquisition") +
  geom_hline(yintercept = 0, linetype =2, color= "#db4437")

dev.off()

#----- External Threats Subfigure -----#

#-- getting the predicted values
external <- ggpredict(m1.2, terms = "lag.roll_avg_hostlev")

#-- plotting

tiff("Fig2External.tiff", units="in", width=5, height=5, res=300)

ggplot(external, aes(x, predicted)) + 
  geom_line(color="#259797") +
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high), fill="#259797", alpha=0.15) +
  theme_bw() +
  ggtitle("External Threats") +
  theme(plot.title = element_text(hjust = 0.5, family="Times"), axis.title = element_text(family="Times"), axis.text = element_text(family="Times")) +
  labs(x= "Average MID Hostility Level", y = "Predicted Probability of Smart Bomb Acquisition") +
  geom_hline(yintercept = 0, linetype =2, color= "#db4437")

dev.off()

#--- Creating Figure 5---#

datainteraction <- data %>%
  rename(Democracy="lagdemoc")

#-- Create special figure theme

theme_special <- function(base_size = 11, base_family = "",
                          base_line_size = base_size / 22,
                          base_rect_size = base_size / 22) {
  # Starts with theme_grey and then modify some parts
  theme_grey(
    #base_size = base_size,
    #base_family = base_family,
    #base_line_size = base_line_size,
    #base_rect_size = base_rect_size
  ) %+replace%
    theme(
      # white background and dark border
      panel.background = element_rect(fill = "white", colour = NA),
      panel.border     = element_rect(fill = NA, colour = "grey20"),
      # make gridlines dark, same contrast with white as in theme_grey
      panel.grid = element_line(colour = "grey92"),
      panel.grid.minor = element_line(size = rel(0.5)),
      # contour strips to match panel contour
      strip.background = element_rect(fill = "white", colour = "grey20"),
      # match legend key to background
      legend.background = element_blank(),
      legend.box.background = element_blank(),
      legend.key = element_blank(),
      complete = TRUE
    )
}

set_theme(
  base = theme_special(), 
  theme.font = "Times", legend.title.face = "plain", legend.item.backcol = "white"
)


#-- Model generation
internal3wayfixedlm <- lm(log_pgm_count_stata~lag.roll_internal_conflict_intensity+lag.treaties_signed+Democracy+lag.roll_avg_hostlev+t+t2+t3+
                            lag.treaties_signed*lag.roll_internal_conflict_intensity*Democracy+factor(ccode) - 1, data=datainteraction)

external3wayfixedlm <- lm(log_pgm_count_stata~lag.roll_internal_conflict_intensity+lag.treaties_signed+Democracy+lag.roll_avg_hostlev+t+t2+t3+
                            lag.treaties_signed*lag.roll_avg_hostlev*Democracy+factor(ccode) - 1, data=datainteraction)

#---- Internal Subfigure -----#

#-- plotting

tiff("Internal3WayInteraction.tiff", units="in", width=5, height=5, res=300)

plot_model(internal3wayfixedlm, type = "pred", terms = c("lag.treaties_signed","lag.roll_internal_conflict_intensity[0,0.1669,2]","Democracy"), mdrt.values = "meansd",
           colors = c("#259797","#208383","#165a5a"), legend.title = "UCDP Internal Conflict Intensity") +
  ggtitle("Internal Threats") +
  labs(x= "ICRC Treaties", y = "Logged Smart Bomb Count (Predicted)") +
  theme(plot.title = element_text(hjust = 0.5), legend.position = "bottom")

dev.off()

#----- External Subfigure -----#

#-- plotting

tiff("External3WayInteraction.tiff", units="in", width=5, height=5, res=300)

plot_model(external3wayfixedlm , type = "pred", terms = c("lag.treaties_signed","lag.roll_avg_hostlev[0,1.156,5]","Democracy"), mdrt.values = "meansd",
           colors = c("#259797","#208383","#165a5a"), legend.title = "Average MID Hostility Level") +
  ggtitle("External Threats") +
  labs(x= "ICRC Treaties", y = "Logged Smart Bomb Count (Predicted)") +
  theme(plot.title = element_text(hjust = 0.5), legend.position = "bottom") 

dev.off()

#--- Creating Figure 6 ---#

#-- read in data
GDPStata <- read.csv('GDPNoTime.csv')

#-- plotting

tiff("GDPLinearPrediction.tiff", units="in", width=5, height=5, res=300)

ggplot(GDPStata, aes(GDPpercapita, Margin)) + 
  geom_line(color="#259797") +
  geom_ribbon(aes(ymin=X.95..Conf., ymax=Interval.), fill="#259797", alpha=0.15) +
  theme_bw() +
  theme(axis.title = element_text(family="Times"), axis.text = element_text(family="Times")) +
  labs(x= "Lagged Logged GDP Per Capita", y = "Linear Prediction of Smart Bomb Acquisition") +
  geom_hline(yintercept = 0, linetype =2, color= "#db4437")

dev.off()

#---Creating Figure 7---#


#----- High Tech Exports Subfigure -----#

#-- read in data
HTEStata <- read.csv('HighTechExportsNoTime.csv')

#-- plotting
tiff("HTEStataPlot.tiff", units="in", width=5, height=5, res=300)

ggplot(HTEStata, aes(HighTechExportsValue, Margin)) + 
  geom_line(color="#259797") +
  geom_ribbon(aes(ymin=X.95..Conf., ymax=Interval.), fill="#259797", alpha=0.15) +
  theme_bw() +
  theme(axis.title = element_text(family="Times"), axis.text = element_text(family="Times")) +
  labs(x= "Lagged Logged High Tech Exports", y = "Linear Prediction of Smart Bomb Acquisition") +
  geom_hline(yintercept = 0, linetype =2, color= "#db4437")

dev.off()

#----- Science and Technology Articles Subfigure -----#

#-- read in data
STAStata <- read.csv('SciTechArticlesNoTime.csv')

#-- plotting

tiff("STAStataPlot.tiff", units="in", width=5, height=5, res=300)

ggplot(STAStata, aes(SciTechArticlesValue, Margin)) + 
  geom_line(color="#259797") +
  geom_ribbon(aes(ymin=X.95..Conf., ymax=Interval.), fill="#259797", alpha=0.15) +
  theme_bw() +
  theme(axis.title = element_text(family="Times"), axis.text = element_text(family="Times")) +
  labs(x= "Lagged Logged Science and Technology Articles", y = "Linear Prediction of Smart Bomb Acquisition") +
  geom_hline(yintercept = 0, linetype =2, color= "#db4437")

dev.off()

#----- Patents Subfigure -----#
#-- read in data
PStata <- read.csv('PatentsNoTime.csv')

#-- plotting

tiff("PStataPlot.tiff", units="in", width=5, height=5, res=300)

ggplot(PStata, aes(PatentsValue, Margin)) + 
  geom_line(color="#259797") +
  geom_ribbon(aes(ymin=X.95..Conf., ymax=Interval.), fill="#259797", alpha=0.15) +
  theme_bw() +
  theme(axis.title = element_text(family="Times"), axis.text = element_text(family="Times")) +
  labs(x= "Lagged Logged Patent Applications", y = "Linear Prediction of Smart Bomb Acquisition") +
  geom_hline(yintercept = 0, linetype =2, color= "#db4437")

dev.off()

#----- IDI Subfigure -----#
#-- read in data
IDIStata <- read.csv('IDINoTime.csv')

#-- plotting

tiff("IDIStataPlot.tiff", units="in", width=5, height=5, res=300)

ggplot(IDIStata, aes(IDIValue, Margin)) + 
  geom_line(color="#259797") +
  geom_ribbon(aes(ymin=X.95..Conf., ymax=Interval.), fill="#259797", alpha=0.15) +
  theme_bw() +
  theme(axis.title = element_text(family="Times"), axis.text = element_text(family="Times")) +
  labs(x= "Lagged Logged IDI", y = "Linear Prediction of Smart Bomb Acquisition") +
  geom_hline(yintercept = 0, linetype =2, color= "#db4437")

dev.off()

#----- Creating Table A.1: Summary Statistics for Robustness Variables -----#


#-- read in data

forrobustness <- data %>%
  select(lag.num_rivals,lag.total_borders,lag.roll_mid_yearly_count,lag.interstate_conflict_intensity,lag.total_internal_conflicts,lag.cinc_5,lag.milex_5,lag.milex_s_i,lag.membership, ccode, year)

additionalvariables <- read.csv('AdditionalVariables1.csv')

forrobustness <- forrobustness %>%
  left_join(additionalvariables)

forrobustness <- forrobustness %>%
  select(-year, -ccode)

#-- create table
stargazer(forrobustness, style = "ajps", summary.stat = c("n", "mean", "sd","min","max"), out="robustsummarystats.tex",
          covariate.labels=c("Lagged Total Rivals", "Lagged Total Borders", "Lagged Average MID Yearly Count",
                             "Lagged UCDP Interstate Conflict Intensity","Lagged UCDP Total Internal Conflicts","Lagged CINC Score",
                             "Lagged COW Military Expenditure", "Lagged SIPRI Military Expenditure", "Lagged Total IGO Membership",
                             "Lagged Log Science and Technology Articles","Lagged Log High Technology Exports","Lagged Log Patent Applications","Defense Relationship Exists","Major Power")
)



